close all
clear
clc
%%%%%%%%%%%%%%%%%%%%%
%     Parameter     %
%%%%%%%%%%%%%%%%%%%%%
tic
a=0.2;                      %lowpass filter not used
slideWindow=300;             %delta function forwarding and tailing window
cutoff=3;                  %peak selection cutoff critias
excludePercentage=0;      %fitting excluding percentage
fitWindow=100;
cutoffForce=10;
stepCutoffMagnet=10;
internalTTest=1;            %Single T test
internalTTestAlpha=0.01;
position_overstretching=12.41;
aa_length=0.5;    % the pesistent length of a amino acid
filelist=importdata('G:\Research_Data\talinr1r3_ivvi_drug\filelist.txt');
only_unfolding=0;% only process the data in the unfolding part of spectrum
decomposecycles=0; % Separate data into different forcecycles by addition of a numerical flag (only the force altering region is outpeed
only_refolding=1;
with_smoothing=1;%whether to smooth rawdata
drift_correction=1; %whether to correct linear drift in long time measurements
forcechange_correction=0;%whether force is changing when steps are detected changes fitting parameters.
forcechange_correction_direct=0;
drift_correction_local=1;
first_points_ignore=500;
iSaveX = @(fname,x) save(fname,'x','-ascii','-double','-tabs');
iSaveX1 = @(fname,x) save(fname,'x');

tic
index=1;
parfor currentfile =1:length(filelist.textdata')
    
    name=char(filelist.textdata(currentfile));
    %% Data processing
    steps_data=[];
    cycledata=[];
    filename=sprintf('%s%s',name,'.txt')
    %c=65/(exp(-(13.5-position_overstretching)/0.36)+0.48*exp(-(13.5-position_overstretching)/1.12));
    time_frame=filelist.data(currentfile,:);
    
    %     datastruct=importdata(filename);
    %     dataraw=datastruct.data;
    %     egxT=dataraw(:,1)';
    %     egMagneticT=dataraw(:,2);
    %     egRawDynT=dataraw(:,3);
    
    
    % Open file
    fid3 = fopen(filename); % Open file
    
    dataraw=textscan(fid3,'%f %f %f %f %f %f %f %f %f %f','HeaderLines',2);
    % Close file
    fclose(fid3);
    
    
    %     datastruct=importdata(filename);
    %     dataraw=datastruct.data;
    egxT_raw=cell2mat(dataraw(1));
    egMagneticT_raw=cell2mat(dataraw(2));
    egRawDynT_raw=cell2mat(dataraw(3));
    
    ind = find(egxT_raw>= time_frame(1) & egxT_raw < time_frame(2));
    name=sprintf('%s%s%d%s%d',name,'_',time_frame(1),'_',time_frame(2));
    
    egxT=egxT_raw(ind);
    egMagneticT=egMagneticT_raw(ind);
    egRawDynT=egRawDynT_raw(ind);
    
    
    
    %plot(egxT); % Check equal spacing
    
    %Mingxi data: A. linear drift remover
    %             B. force altering regions
    %             C. points for fitting selection
    
    if drift_correction
        linearFitP=polyfit(egxT(egMagneticT==min(egMagneticT)),egRawDynT(egMagneticT==min(egMagneticT)),1);
        egfittedDriftT=polyval(linearFitP, egxT)-polyval(linearFitP, min(egxT));
        disp('Determining the drift\n');
        name
        linearFitP % display the drift
        egRawDynTShift=egRawDynT-egfittedDriftT;
    else
        egRawDynTShift=egRawDynT;
    end
    
    customDyn=egRawDynTShift;
    
    x=egxT;
    
    BState=2
    
    BStateMin = 2;
    BStateMax = 3+BState;
    %set(handles.BStateNumS,'Min',BStateMin);
    %set(handles.BStateNumS,'Max',BStateMax);
    %set(handles.BStateNumS,'SliderStep',[1/(BStateMax-BStateMin) 1/(BStateMax-BStateMin)]);
    
    %set(handles.BStateNumS,'Value',BState);
    
    %set(handles.BStateNumE,'String',sprintf('%d', round(get(handles.BStateNumS,'Value'))));
    
    %Initialization
    iniTMatrix = ones(BState).*(1/BState);
    iniEPara = [prctile(customDyn,(1:2:2*BState)/(2*BState)*100)' ones(BState,1).*std(customDyn)];
    pFunction = @(x,y)normpdf(x,y(1),y(2));
    
    
    
    TMatrixFix = zeros(BState);
    EParaFix = zeros(BState, 2);
    [TMatrix, EPara, logLikelihood] = MHMMTRAIN_mingxi_mex(customDyn, iniTMatrix, iniEPara, TMatrixFix, EParaFix);
    
    %set(handles.BTMatrixT,'data',TMatrix);
    %set(handles.BEParaT,'data',EPara);
    %set(handles.BTMatrixFixT,'data',logical(TMatrixFix));
    %set(handles.BEParaFixT,'data',logical(EParaFix));
    
    iSaveX(sprintf('%s%s%d%s%d%s',name,'_hmm_TMatrix-',slideWindow,'slide',BState,'state.dat'),TMatrix);
    iSaveX(sprintf('%s%s%d%s%d%s',name,'_hmm_EPara-',slideWindow,'slide',BState,'state.dat'),EPara);
    iSaveX(sprintf('%s%s%d%s%d%s',name,'_hmm_logLikelihood-',slideWindow,'slide',BState,'state.dat'),logLikelihood');
    
    [FilDyn, FilDynState] = MHMMVITERBI_mingxi_mex(customDyn, TMatrix, EPara);
    
   
    iSaveX(sprintf('%s%s%d%s%d%s',name,'hmm_states',slideWindow,'slide',BState,'state.dat'),[x FilDynState FilDyn customDyn]);
    
    %export_fig gcf figure_name -png;
    % Plotting
    %             figure()
    %             plot(handles.FilDynAxes,FilDyn);
    %             axis(handles.FilDynAxes,[1 length(FilDyn) min(FilDyn)-.1*abs(max(FilDyn)-min(FilDyn)) max(FilDyn)+.1*abs(max(FilDyn)-min(FilDyn))]);
    %
    %             plot(handles.BLikelihood,logLikelihood,'b');
    
    % Updating
    %             handles.metricdata.FilDyn = FilDyn;
    %             handles.metricdata.FFLAG = 1;
    
%% calculate the states of the life time assuming two state at the moment  
statedata=[x FilDynState FilDyn customDyn];

diffstates_lifetime=diff(statedata(diff(statedata(:,2))~=0,1:2));
logdist_state1=log((diffstates_lifetime(diffstates_lifetime(:,2)==-1,1)));
logdist_state2=log((diffstates_lifetime(diffstates_lifetime(:,2)==1,1)));

iSaveX(sprintf('%s%s%d%s%d%s',name,'_hmm_loglifetime_state1-',slideWindow,'slide',BState,'state.dat'),logdist_state1);
iSaveX(sprintf('%s%s%d%s%d%s',name,'_hmm_loglifetime_state2-',slideWindow,'slide',BState,'state.dat'),logdist_state2);

[loghist_state1_elements,loghist_state1_center]=hist(logdist_state1,30);
[loghist_state2_elements,loghist_state2_center]=hist(logdist_state2,30);

fitfunction=@(a,x)  sqrt(a(2)*exp(x-a(1)-exp(x-a(1))));
state1model=NonLinearModel.fit(loghist_state1_center,sqrt(loghist_state1_elements),fitfunction,[1 max(sqrt(loghist_state1_elements))]);
state2model=NonLinearModel.fit(loghist_state2_center,sqrt(loghist_state2_elements),fitfunction,[0.5 max(sqrt(loghist_state2_elements))]);

iSaveX1(sprintf('%s%s%d%s%d%s',name,'_hmm_loglifetime_state1_fitting_model-',slideWindow,'slide',BState,'state.mat'),state1model);
iSaveX1(sprintf('%s%s%d%s%d%s',name,'_hmm_loglifetime_state2_fitting_model-',slideWindow,'slide',BState,'state.mat'),state2model);
g=figure;
range1=min(loghist_state1_center):0.01:max(loghist_state1_center);
stairs(loghist_state1_center,sqrt(loghist_state1_elements),'g','LineWidth',2);
hold on;
plot(range1,predict(state1model,range1'),'r','LineWidth',2,'LineSmoothing','on');
figure_name=sprintf('%s%s%d%s%d%s',name,'_hmm_states_loghistogram-',slideWindow,'slide',BState,'state.fig');
hold on;
range2=min(loghist_state2_center):0.01:max(loghist_state2_center);
stairs(loghist_state2_center,sqrt(loghist_state2_elements),'b','LineWidth',2)
hold on;
plot(range2,predict(state2model,range2'),'r','LineWidth',2,'LineSmoothing','on');
legend(sprintf('%s%0.2f%s','state1 ',EPara(1),' nm'),sprintf('%s%0.2f%s','state1 fit tau = ',exp(state1model.Coefficients.Estimate(1)),'s'),sprintf('%s%0.2f%s','state2 ',EPara(2),' nm'),sprintf('%s%0.2f%s','state2 fit tau = ',exp(state2model.Coefficients.Estimate(1)),' s'),'Location','Best');
%figure_name1=sprintf('%s%s%d%s%d%s',name,'_hmm_state2_loghistogram-',slideWindow,'slide',BState,'state.fig');
saveas(g,figure_name,'fig');
    
end

toc

